clear; close all; clc;

%% 1. 模型参数
alpha = 5;
sigma = 0.8;
mu    = 0.1;
epsilon = 0.1;

% 为了充分脱离瞬态，设置较长迭代次数
N_total = 1000;  % 总迭代次数
N_cut   = 100;   % 去除瞬态

% 扫描范围：phi0 in [-4, 6], k in [0, 0.4]
phi_min = -4;  phi_max = 6;  Nphi = 151;
k_min   = 0;   k_max   = 0.4; Nk   = 151;

phi_values = linspace(phi_min, phi_max, Nphi);
k_values   = linspace(k_min, k_max, Nk);

%% 2. 预分配矩阵存储平均尖峰数
avgSpike_map = zeros(Nphi, Nk);

%% 3. 双参数扫描
parfor i_phi = 1:Nphi
    phi0_now = phi_values(i_phi);
    tmp_avgSpike = zeros(1, Nk);
    for i_k = 1:Nk
        k_now = k_values(i_k);
        % 调用 Rulkov 模型，初始状态在 Rulkov.m 中固定为 (0.01, 0, -0.01, 0)
        [X1_save, ~, ~] = Rulkov(alpha, sigma, mu, k_now, phi0_now, N_total, epsilon, N_cut);
        % 取去除瞬态后的有效数据
        X_eff = X1_save((N_cut+1):end);
        
        % 对信号进行平滑处理（采用移动中值滤波，窗口大小为5）
        X_smooth = smoothdata(X_eff, 'movmedian', 5);
        
        % 计算该平滑信号的平均尖峰数
        asn = AvgSpikeNum(X_smooth);
        
        % 如果平均尖峰数大于20，则统一记为20
        if asn > 20
            asn = 20;
        end
        
        tmp_avgSpike(i_k) = asn;
    end
    avgSpike_map(:, i_phi) = tmp_avgSpike;
end

%% 4. 绘制二维平均尖峰数分布图
figure;
imagesc(k_values, phi_values, avgSpike_map);
axis xy;  % 保证 y 轴方向从下到上
xlabel('k');
ylabel('\phi_0');
title('(c) 平均尖峰数');
colormap jet;
cb = colorbar;
cb.Label.String = '平均尖峰数';
ylim([phi_min, phi_max]);  % 保证 y 轴显示范围为 [-4, 6]
clim([1, 20]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function asn = AvgSpikeNum(X)
% 计算时间序列 X 中每个 burst 的平均尖峰数
%
%   1. 利用 findpeaks 检测所有局部峰值；
%   2. 根据相邻峰值的索引差分，设定阈值（这里取 1.8 倍的平均间隔）区分 burst；
%   3. 对每个 burst 统计尖峰个数，求平均值。
%
% 若检测不到 burst，则返回 1。

    % 1. 检测所有峰值及其索引
    [~, locs] = findpeaks(X);
    if isempty(locs) || numel(locs) < 2
        asn = 1;
        return;
    end
    
    % 2. 计算相邻峰值间隔和平均间隔
    diffs = diff(locs);
    avgDiff = mean(diffs);
    % 设置阈值：相邻峰间隔大于 1.8*avgDiff 认为是不同 burst 的分界
    threshold = 1.8 * avgDiff;
    
    bursts = {};  % 存储每个 burst 的峰值索引
    currentBurst = locs(1);
    burstCount = 1;
    
    for i = 2:length(locs)
        if diffs(i-1) > threshold
            bursts{burstCount} = currentBurst; %#ok<AGROW>
            burstCount = burstCount + 1;
            currentBurst = locs(i);
        else
            currentBurst = [currentBurst, locs(i)]; %#ok<AGROW>
        end
    end
    bursts{burstCount} = currentBurst;
    
    % 3. 计算每个 burst 的峰值数，并求平均
    counts = cellfun(@length, bursts);
    asn = mean(counts);
    
    % 最少记为 1
    if asn < 1
        asn = 1;
    end
end
